pacman::p_load(tidyverse,data.table,broom,raster,latex2exp,DescTools,
               rcartocolor,ggpubr,ggpattern)
rm(list=ls())
################################################################################


##################### Emissions footprints by commodity 

#By stage of the supply chain
df = read_csv('../../data/environmental/raw/GHG-emissions-by-life-cycle-stage-OurWorldinData-upload.csv',show_col_types = FALSE)
setnames(df, old=c('Food product','Land use change','Animal Feed'), new=c('product','LUC','Feed'))
animals = c('Beef (beef herd)', 'Fish (farmed)', 'Poultry Meat','Pig Meat','Eggs')
plants = c('Rice', 'Tofu', 'Potatoes','Wheat & Rye (Bread)','Other Pulses','Maize (Meal)')
df <- df %>% filter(product  %in% animals | product  %in% plants)
df <- df %>% mutate(Postfarm = Processing + Transport +  Packging + Retail )
df <- df[, names(df) %in% c("product","LUC","Feed","Farm","Postfarm")]
df$product <- reorder(df$product, rowSums(df[-1]))
df <-gather(df, stage, emissions, LUC:Postfarm, factor_key=TRUE)
df=setDT(df)
df[stage=='LUC',]$stage = "Land use change"
df[stage=='Postfarm',]$stage = "Post-farm"
df[stage=='Farm',]$stage = "On-farm"
df[stage=='Feed',]$stage = "Feed"
df$stage <- factor(df$stage, levels = c("Post-farm",
                                        "On-farm",
                                        "Feed",
                                        "Land use change"))

figure1a = ggplot(df, aes(x = product, y=emissions, fill=stage, pattern=stage)) + 
  geom_bar_pattern(stat='identity', linewidth=0.1, colour='gray30', 
                   pattern_color= "gray30",
                   pattern_density = 0.0025,
                   pattern_spacing = 0.04,
                   pattern_key_scale_factor = 0.5) +
  labs(y=TeX(r'(kg of CO$_2$e per kg of product)'), title='A. Emissions footprints and sources' ) +
  theme_minimal() +
  theme(plot.title=element_text(size=9, hjust=1, vjust=1),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size=8),
        axis.text.x = element_text(size=7),
        axis.text.y = element_text(size=7),
        legend.title = element_blank(), 
        legend.position='inside', 
        legend.position.inside = c(0.7, 0.3),  
        legend.key.size = unit(0.25, "cm"),
        legend.direction = "vertical", legend.text = element_text(size=7),
        aspect.ratio = 0.7,
        panel.grid.minor = element_blank()
  ) +
  scale_x_discrete(labels=c("Beef (beef herd)"="Beef",
                            "Pig Meat"="Pork",
                            "Poultry Meat"="Chicken",
                            "Other Pulses"="Legumes",
                            "Tofu" = "Tofu (soybean)",
                            "Maize (Meal)" = "Maize",
                            "Wheat & Rye (Bread)" = "Wheat & Rye" )) +
  scale_pattern_manual(values=c('wave','stripe','circle','none')) +
  scale_fill_manual(values=c("#999999","#D55E00","#E69F00","#009E73")) +
  coord_flip()
figure1a_bw <- figure1a + scale_fill_grey()


##################### Carbon maps

raster_file = raster(paste0("../../data/environmental/raw/aboveground.asc"))
raster_file_pts <- rasterToPoints(raster_file, spatial = TRUE)
raster_file_df  <- data.frame(raster_file_pts)
rm(raster_file_pts, raster_file)
raster_file_df <- raster_file_df %>% filter(aboveground > 0)
raster_file_df$aboveground <- Winsorize(raster_file_df$aboveground, maxval = 750)
figure1b = ggplot() +
  geom_tile(data = raster_file_df , aes(x = x, y = y, fill = aboveground)) + 
  scale_fill_carto_c(name = 'tC/ha', palette = "Mint",  na.value = "transparent",  guide='colourbar') +
  coord_fixed() +
  theme_void() + 
  labs(fill = 'tC/ha', title='B. Carbon intensity of land') + 
  theme(legend.position = "right",  legend.direction = "vertical", 
        plot.title=element_text(size=9, hjust=-0.5, vjust=1),
        legend.title = element_text(size=8),
        legend.text = element_text(size=7),
        legend.key.size = unit(0.6, "cm"),
        legend.key.width = unit(0.25,"cm"))

figure1b_bw = ggplot() +
  geom_tile(data = raster_file_df , aes(x = x, y = y, fill = aboveground)) + 
  scale_fill_gradient(low = "gray90", high = "gray20") +
  coord_fixed() +
  theme_void() + 
  labs(fill = 'tC/ha', title='B. Carbon intensity of land') + 
  theme(legend.position = "right",  legend.direction = "vertical", 
        plot.title=element_text(size=9, hjust=-0.5, vjust=1),
        legend.title = element_text(size=8),
        legend.text = element_text(size=7),
        legend.key.size = unit(0.6, "cm"),
        legend.key.width = unit(0.25,"cm"))

#Combined figure
ggarrange(figure1a,figure1b, nrow=1,ncol=2)
ggsave("../../output/figures/figure1.pdf", width = 5.5, height = 2, units='in')
ggsave("../../output/figures/figure1.eps", width = 5.5, height = 2, units='in')

ggarrange(figure1a_bw,figure1b_bw, nrow=1,ncol=2)
ggsave("../../output/figures/figure1_bw.pdf", width = 5.5, height = 2, units='in')
ggsave("../../output/figures/figure1_bw.eps", width = 5.5, height = 2, units='in')


##################### Appendix figures - emissions footprints by commodity 

#Per 1000 kcal
df = read_csv('../../data/environmental/raw/ghg-kcal-poore.csv',show_col_types = FALSE)
setnames(df, old=c('GHG emissions per 1000kcal (Poore & Nemecek, 2018)'), new=c('emissions'))
animals = c('Beef (beef herd)', 'Fish (farmed)', 'Poultry Meat','Pig Meat','Eggs')
plants = c('Rice', 'Tofu (soybeans)', 'Potatoes','Wheat & Rye','Other Pulses','Maize')
df <- df %>% filter(Entity  %in% animals | Entity  %in% plants )
df <- df %>% mutate(type = ifelse(Entity  %in% animals, "Animal based", "Plant based"))
figure_kcal <- ggplot(df, aes(reorder(Entity,emissions), emissions, fill=type, pattern=type)) + 
  geom_bar_pattern(stat='identity', linewidth=0.35, colour='black', 
                   pattern_color= "gray30",
                   pattern_density = 0.0025,
                   pattern_spacing = 0.03,
                   pattern_key_scale_factor = 0.75) +
  labs(y=TeX(r'(kg of CO$_2$e per 1000 kcal)')) +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(margin = margin(t = 10), size=14),
        legend.title = element_blank(), legend.position='inside', legend.position.inside = c(0.825, 0.15),  legend.direction = "vertical", legend.text = element_text(size=14),
        aspect.ratio = 1,
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(size=17),
        axis.text.y = element_text(size=17)) +
  scale_pattern_manual(values=c('none','stripe')) +
  scale_fill_manual(values=c("#D55E00","#009E73")) +
  scale_x_discrete(labels=c("Beef (beef herd)"="Beef","Pig Meat"="Pork",
                            "Poultry Meat"="Chicken","Other Pulses"="Legumes")) +
  coord_flip()

#Per kg of protein
df = read_csv('../../data/environmental/raw/ghg-per-protein-poore.csv',show_col_types = FALSE)
setnames(df, old=c('GHG emissions per 100g protein (Poore & Nemecek, 2018)'), new=c('emissions'))
animals = c('Beef (beef herd)', 'Fish (farmed)', 'Poultry Meat','Pig Meat','Eggs')
plants = c('Rice', 'Tofu (soybeans)', 'Potatoes','Wheat & Rye','Other Pulses','Maize')
df <- df %>% filter(Entity  %in% animals | Entity  %in% plants )
df <- df %>% mutate(type = ifelse(Entity  %in% animals, "Animal based", "Plant based"))
df$emissions = df$emissions*10 
figure_protein <- ggplot(df, aes(reorder(Entity,emissions), emissions, fill=type, pattern=type)) + 
  geom_bar_pattern(stat='identity', linewidth=0.35, colour='black', 
                   pattern_color= "gray30",
                   pattern_density = 0.0025,
                   pattern_spacing = 0.03,
                   pattern_key_scale_factor = 0.75) +
  labs(y=TeX(r'(kg of CO$_2$e per kg of protein)') ) +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(margin = margin(t = 10), size=14),
        legend.title = element_blank(), legend.position='inside', legend.position.inside = c(0.825, 0.15),  legend.direction = "vertical", legend.text = element_text(size=14),
        aspect.ratio = 1,
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(size=17),
        axis.text.y = element_text(size=17) ) +
  scale_pattern_manual(values=c('none','stripe')) +
  scale_fill_manual(values=c("#D55E00","#009E73")) +
  scale_x_discrete(labels=c("Beef (beef herd)"="Beef","Pig Meat"="Pork",
                            "Poultry Meat"="Chicken","Other Pulses"="Legumes")) +
  coord_flip()

ggarrange(figure_kcal,figure_protein,nrow=1,ncol=2)
ggsave(paste0("../../output/figures/appendix/figure10.pdf"), width = 12.5, height = 5)
ggsave(paste0("../../output/figures/appendix/figure10.eps"), width = 12.5, height = 5)

